Mixing-induced activity in open flows 
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Abstract. We develop a theory describing how a convectively unstable active held 
in an open flow is transformed into absolutely unstable by local mixing. Presenting the 
^yy ' mixing region as one with a locally enhanced effective diffusion allows us to find the 

p I . linear transition point to an unstable global mode analytically. We derive the critical 

£H \ exponent that characterizes weakly nonlinear regimes beyond the instability threshold 

and compare it with numerical simulations of a full two-dimensional flow problem. The 
obtained scaling law turns out to be universal as it depends neither on geometry nor 
on the nature of the mixing region. 

> ' 

^ I PACS numbers: 47.54.-r; 47.70.-n; 89.75.Kd 

m _ 

! 1. Introduction 

In many geophysical and laboratory flows active chemical and biological processes occur. 
These active processes are crucially dependent on the nature of the flow. Especially 
important from the theoretical point of view becomes understanding the role of mixing. 
As we demonstrate below, even localized region of strong mixing introduced in the 
laminar flow is able to significantly change the overall picture of activity. Because the 
active processes accompanied by mixing are widespread, the outcomes of our analysis 
are highly relevant to practically important problems varying from chemical reactions 
in micro-mixers to plankton growth in ocean, see [1] and references therein. 

What turns out to be generic, the active processes quite often occur in an open 
rather than in a closed geometry. Here the main issue is whether the throughflow is 
stronger or weaker than the activity. One has to compare the velocity of the throughflow 
with the velocity of the activity spreading due to diffusion. If the throughflow is stronger, 
the activity is blown away like a candle flame in a strong wind, in the opposite case 
a sustained activity can be observed [21 El H]. This simple picture is valid, however, 
only for homogeneous media. Often additional vortexes are superimposed on a constant 
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throughflow, due to, e.g., wakes behind islands in ocean currents or mixing enforced 
by revolving fan blades in laboratory experiments. We want to study under which 
conditions such an additional kinematic mixing in a strong open flow can lead to a 
transition to a sustained activity, and to characterize this transition quantitatively. 

Our main model is a reaction-advection-diffusion equation for the dimensionless 
concentration of an active scalar field u(r, t) 

— + [V + W(r, t)} ■ Vw = A) V 2 m + au(l - u p ) . (1) 

Here V = (V, 0, 0) is a constant throughflow in x-direction, Do is a molecular diffusion 
of the scalar field. Activity is assumed to be of the simplest form: a linear growth with 
rate a with a saturation at u — 1. The nonlinearity index p is typically integer (1 or 2) 
for chemical reactions, while for biological populations a wide range of values of p has 
been recently reported [5]. Mixing is described by a spatially localized incompressible 
velocity field W(r, t), its intensity is denoted as W. Note that in the absence of fluid 
flow Eq. (pp) is reduced to the famous Kolmogorov-Petrovsky-Piskunov-Fisher (KPPF) 
model of an active medium with diffusion (see, e.g., [6] for original references, analysis, 
and applications of KPPF), while for a = Eq. (CQ) describes a linear evolution of a 
passive scalar in a flow. Model can be used for the description of biological activity, 
where u is, e.g., concentration of a growing plankton, advected by oceanic currents 
[7]; for a possible laboratory realization see recent experiments with an autocatalytic 
reaction in a Hele-Shaw cell with a throughflow [8]. 

In the absence of flow, the diffusion causes the active state to spread forming 
eventually a front with velocity Vf = 2y/aD |9J. Thus, for vanishing mixing W = 0, 
the activity is blown away provided V > Vf. For this parameter range the instability 
in Eq. ([I]) is convective and in the absence of external sources, no activity is observed 
in the medium. A nontrivial state is, however, possible if there is a localized source of 
the field u: then a growing tail stretches from this source in the downstream direction, 
where it eventually saturates at u — 1. The linearized problem with a point (5-function) 
source of intensity e can be readily solved, yielding 

u(x) = e(V)- 1 exp[xV/2D \ exp(-\x\V/2D Q ) (2) 

in one-dimensional setups and 

u(r) = e(2nD )- 1 exp[xV/2D ]K (\r\V/2D ) (3) 

in two dimensions (where K is the modified Bessel function, V = (V 2 — Vj) 1 ^ 2 ). Note 
that in both solutions u ~ exp(/i T x) as x — > ±oo, where p,± = (2D ) ^ 1 (V ± V). 

In this paper we demonstrate, that beyond some critical intensity W cr , a localized 
mixing W(r, t) turns the convective instability locally into the absolute one, which 
results in a stationary (in statistical sense) profile of u (see Fig. [1] for a sketch of the 
profile and Figs.H]and[7|below for numerical examples). Beyond criticality W > W cr , the 
mixing region acts as an effective source of the field, in Fig. [1] this region is denoted as a 
"source." A "tail" where the field grows exponentially as in (J2J), ([3]) extends downstream 
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of the source. Further downstream stretches the "plateau" domain where u = 1. Our 
main quantitative result, obtained by matching solutions in these three domains, is the 
critical exponent j3 relating the intensity of the effective source e e ff with the mixing 
intensity: e ef f ~ (W - W CT Y . 

To develop the theory we use the concept of global modes pHl [HI [12]. In this 
concept a self-sustained non-advected pattern arises due to inhomogeneities of the 
system. Typically, one considers inhomogeneities of the growth rate a: if a = a(r) has 
a hump where locally the front velocity is large Vj oc > V, then a global mode appears, 
located at this hump and downstream. In this paper we are interested in another, 
mixing-based mechanism of a global mode appearance. It can be easily understood 
if the concept of effective diffusion (see, e.g., [13J) is used to describe the mixing 
term in (CQ). In this approach, we phenomenologically introduce effective diffusivity 
D(r) = Dq + D m i X (r) that accounts for the coarse grained mixing dynamics, and write 
instead of §B) an equation with a non-homogeneous diffusion 
3u 

— + V • Vu = V[£>(r)Vw] + au(l - u p ) . (4) 

C/ v 

A hump of diffusivity D(r) leads to an increase of the local front velocity V/, and one 
expects that when the front propagation prevails over the throughflow, a stationary 
global mode can appear, producing a mixing-induced sustained structure. We focus 
on a geometry shown in Fig. [T](a): in a constant open flow there is a localized region 
of strong mixing, which, as we will see below, is not necessarily chaotic or turbulent. 
The theory below will be developed for a one-dimensional case, which is relevant, e.g., 
for flows in a micropipe; the results will be supported by numerical simulations of two- 
dimensional flows. We restrict ourselves to this case because of computational simplicity, 
and also because two-dimensional flows are relevant for many geophysical and laboratory 
(especially in microfluidics) experimental situations. 



2. Linear stability analysis 

2.1. Effective diffusion model 

We start with a linear analysis of a one-dimensional situation, described by the linearized 
at u = Eq. gj): 



® u _|_ y® u ® 
dt dx dx 



+ au . 



(5) 



We look for an exponentially growing in time solution and with an ansatz u(x,t) 
exp[At + f x z(£)d£\ obtain 

dD(x) 



dz 
dx 



V 



dx 



a — A 



■z — 



(6) 



As \x\ 



D(x) D(x) 

oo we have a homogeneous medium with D = Dq, here the solution should 



tend to values z± = (2D ) 1 (V ±^V 2 — 4(a — X)D ) at which the r.h.s. of §6§ vanishes. 
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(a) "source" "tail" "plateau" 




Figure 1. (a): Quasi-one-dimensional open flow with a localized mixing zone. 
Panels (b) and (c) illustrate the construction of the nonlinear global mode, they show 
qualitative profiles u(x) and z(x) = in the linear approximation at the criticality 

(dashed line) and in nonlinear regime slightly beyond criticality (full line). In the 
latter case the profile is nearly linear for x < x r but deviates due to nonlinear terms 
for x > x r , see discussion of Eq. (Tl6)) . Regions "source," "tail," and "plateau" are 
explained in the text. 



More precisely, as x — > — oo we have z — > z+ and as x — > oo we have z — > z°_ . With these 
two boundary conditions one easily finds the solution of ([6]) numerically, matching at 
z = integrations starting at large \x\ from the values z±. In a particular analytically 
solvable case of a piecewise-constant diffusivity: D = Do for \x\ > I and D = D\ > Dq 
for \x\ < I, one can perform the integration analytically and obtain the equation for the 
growth rate A: 



V 2 - 


-4(a- 


A)A) 


4(a 


-A)£>! 


-V 2 



I — — / 1 arctan \ I -J- -\— — ^— ^ . (7) 

y/A(a - X)D 1 - V 2 1 " x "° 1 

The value A = corresponds to the onset of global instability, in this case (J7J) gives the 
relation between the critical values l cr and -D lcr : 



2D la . V^-Vf 

= vm^w arctan V 4« Dur - ■ (8) 

From (jSJ) it follows that Di cr — > oo as l cr — > l m in — V /{2a). In other words, there exists 
a minimal size of the mixing region, so that for I < l m in even a very strong mixing, 
with a very large effective diffusion, cannot create a global mode (the same is true for a 
smooth profile of D(x); note also that the size of the mixing region is not limited from 
above). This is in contrast to the situation when the global mode is induced by a local 
hump of the growth rate a (cf. [15]): here one can obtain instability even when a(x) is 
highly localized (a delta-function), a global mode then looks as in (T5]), ([3]). 
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A similar analysis can be performed for a two-dimensional inhomogeneous linear 
reaction-advection-diffusion equation 

^g + V— = V[D{r)Vu} + au ) (9) 

where we assume D = D\ for r < R and D = Dq for r > R. Presenting the solutions 
in the inner and outer domains as Ui = ipi(r, 6) exp[(2.Di) _1 W cos 6* + Xt] for r < R 
and u = ip (r, 0) exp[(2_D ) ~ x Vr cos 9 + At] for r > R, we obtain equations for ip ij0 
whose solutions can be written down as series in Bessel functions J m and K m . Matching 
these series at r = R leads to rather cumbersome expressions in terms of expansion 
coefficients. As a result, the eigenvalue problem reduces to a matrix equation which we 
solved numerically. In Fig. [2] we show the critical line (corresponding to the condition 
A = 0) on the plane (R,Di) for V/Vf = 3/2. One can see, that similar to the one- 
dimensional case, there exists a minimal radius of the higher diffusivity spot, so that for 
R < Rmin the global mode cannot arise. 

1000 
o 100 
^ 10 



1 

1 2 3 4 5 10 20 

Ry/a/Do 

Figure 2. The critical line for model ((9]). Dashed line denotes the minimal radius for 
the high-diffusive spot yUi/DoR m in ~ 1-905. 




2.2. Non-uniform flow 

Now we discuss the linear stability not in the framework of the effective diffusion 
model (jl]), but in the full reaction-advection-diffusion problem as in Eq. ([!]). After 
the linearization we arrive at a linear stability problem 

fin i 

— +[\ + W(v,t)]-Vu = D V 2 u + au, (10) 

which is non-stationary if the velocity field W is time dependent. For a time- independent 
W the stability is defined by the growth rate A as before. For time dependent fields W, 
the proper way to determine the stability is to calculate the largest Lyapunov exponent 
(LE) A = (| In 1 1«||). This can be done numerically, as described in [16]. Noteworthy, 
in this consideration we are not restricted to a deterministic flow, as the LE can be 
calculated also for a randomly or chaotically time-dependent field W(r, t). 
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We first calculate the growth rate A for a linearized quasi-one- dimensional reaction- 
advection-diffusion equation ([TO]) subject to a constant open flow with a superimposed 
stationary vortex, described by a stream function 

^(x, y) = Vy + W[cos{27ry) - 1] exp(-x 2 6" 2 ) . (11) 

We fix parameters V = 1, b = 1 and evaluate the growth rate for different molecular 
diffusion constants Dq and vortex intensities W (Fig. [3]). Note that the parameter a 
simply shifts A, therefore we plot A — Ao, where Xq = a — V 2 /(AD ) is the growth rate 
for a non-mixed flow with W = 0. We see that the mixing- induced enhancement of field 
growth is mostly pronounced for small diffusion and saturates at W ~ 0.5. 




Figure 3. Growth rate characterizing stability of the global mode mixed by the flow 
with a stationary vortex (1 1 lj) . 

The role of a mixing vortex can be understood in the following way. In the case of 
a passive advection-diffusion process and for a relatively small molecular diffusivity, the 
field is trapped by the vortex. In the limit of vanishing diffusivity the field u outside 
the vortex is blown away, whereas everything inside the vortex is trapped and cannot 
escape for an infinitely long time, as e.g. in [17]. As a result, a pattern in the form of a 
cloud arises. Because of diffusion, such pattern is no longer stable and its concentration 
gradually decays. Hence, under an advection-diffusion process, the cloud exists only for 
a limited time. However, if we now take into consideration activity, this time is spent 
by the active field to grow, which is enough to compensate the loss caused by diffusion. 
Thus, a self-sustained pattern is born, Fig. HI 

Next we calculate the LE for a linearized two-dimensional reaction-advection- 
diffusion equation (TTO!) with a constant open flow and a superimposed oscillating vortex, 
described by the stream function 

ty 2 (x,y,t) = Vy + Wexp{-(x 2 + y 2 )R~ 2 }cos(ujt) . (12) 

Here we fix V — 1, R — 1 and u = 2 and calculate the LE A, see Fig. [51 As before, 
we plot A — Ao, where Aq = a — V 2 / (4/Aj) is the LE for a non-mixed flow with W = 0. 
Again, the mixing-induced enhancement of field growth is mostly pronounced for small 
diffusion. In contrast to the previous case the dependence on W is non-monotonic and 
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Figure 4. An example of the global mode in the quasi-one-dimensional flow (fTTj) with 
V = 1, W = 0.45, evaluated for 6 = 1, -Do = 0.01, a — 1.5, p = 1. Higher and lower 
values of the active field concentration are plotted with the lighter and darker colours, 
respectively. 




Figure 5. Lyapunov exponent characterizing stability of the global mode mixed by 
periodically blinking vortex (fl~2|) . 



has a maximum at W ~ 3. This is the mixing strength at which a chaotic saddle [TJ in 
the Lagrangean particle trajectories appears. A further increase of the vortex intensity 
does not lead, however, to a significant growth of the LE. 

3. Nonlinear analysis. A universal scaling law 

Now we develop a nonlinear theory of the global mode. It is clear that the nonlinear 
saturation stops the exponential growth of a slightly supercritical linear mode and leads 
to a nonlinear solution with a finite amplitude. Our aim is to describe the dependence 
of this amplitude on the deviation from criticality. First we notice, that the very notion 
of the amplitude is here nontrivial. Indeed, the nonlinear solution looks as in Fig. [l](b) 
(cf. Fig. [7] below); it saturates to u = 1 in the downstream direction. However, outside 
the mixing region the field looks like a solution caused by a localized field source. Thus, 
we can take the effective intensity of this source e e // 3 which is proportional to the 
characteristic field amplitude in the mixing region w(0) [see relations (j2J), ©], as the 
order parameter of the transition. The deviation from the criticality we will measure 
with the growth rate A, for which holds A oc W — W cr in the full model flTJ or A oc D — D cr 
in model (j4j). 

We will consider the simplest possible setup, namely the nonlinear modification of 



Mixing-induced activity in open flows 



8 



one-dimensional Eq. (151) : 

+ au{l-u p ). (13) 



<9t cte <9x 



We look for a stationary global mode u(x), and rewrite this equation as the system 

_ = _; 2+ _^ 2 __ + _ u ,, 1 (14) 

— = . (15) 

ax 

We consider this system separately in two spatial domains. The first, linear region, 
includes the inflow and the mixing domains ("source" in Fig. [1]): — oo < x < x r , where 
the field u(x) remains small. In the second, outflow region x r < x < oo, the field 
u further grows ("tail") and nonlinearly saturates ("plateau"). In the linear region, 
because of smallness of the field, we can neglect u v in ( |T4|) . thus we obtain an equation 
similar to ([6]). The only difference is that because we look for a stationary solution, in 
(114j) the term ~ A is absent. Near the criticality, where A is small, we can consider this 
term as a perturbation, therefore the solution of (|14p in the linear region is close to the 
solution of Eq. (jSJ); it has the asymptotic z — > « + as x — > — oo. Due to the perturbation 
term oc A, at the right border of the linear region z deviates from /i_: the deviation 
/U_ — z(x r ) is proportional to A, and, thus, to D — D cr . At x r the field u is small and 
u(x r ) oc u(0). 

Next we consider full equations (JHJ), (|15p in the nonlinear region x > x r . Here 
the solution should tend as x — > oo to the saddle fixed point u — 1, z = 0. Thus, 
starting integration from large values of x in the negative direction, we have to follow 
the stable manifold of this saddle and match this solution at x = x r with the obtained 
above. Because the value to be matched z(x r ) is very close to in the region where 
the solution (z,u) approaches (z(x r ),u(x r )) we can write z 2 ~ — 2{i_Az to obtain 

-^-Az = U + - fi_)Az - ^-vP(x r )e ptl -^- x ^ . (16) 
ax D 

Here, since Az = /i_ — z(x) is small, we have approximated the solution of (TT5T) as 

u ~ u(x r )e^-( x ~ Xr \ Because linear inhomogeneous Eq. ([TBI) is solved in the negative 

in x direction, the solution follows the slowest exponent: Az oc exp[7(x — x r )], where 

7 = min( / u + — 

At the criticality, the region of validity of the exponential solution Az oc exp[7(x — 
x r )] becomes very large. Thus it is dominant for small deviations from criticality A, 
therefore we can estimate the coordinate x s at which the field u saturates (i.e., we reach 
the state u ~ 1 and 2^0) from the relations above: from — //_ rs (z(x r ) — n^)e l( - Xs ~ Xr ^ it 
follows {x s — x r ) ~ — 7 _1 ln(/i_ — ^(x r )). Substituting this in the expression for u(x), we 
obtain u(x r ) oc (//_ — z(x r )) M ^/ 7 . Now we take into account that — z(x r ) oc D — D cr , 
and, because the evolution of u in the interval < x < x r only weakly depends on the 
criticality, e e ff = u(0) oc u(x r ). The final expression for the scaling law of the amplitude 
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Figure 6. The critical exponent calculated for the model |T]) (symbols), compared 
with theoretical prediction (fT8|) (lines). One can clearly see the crossover between two 
regimes of the field saturation in dependence on parameter a, the latter is related to v 
in |18|) via v ~ a~ 1/2 . 



of the global mode thus reads 



E ef f OC , 



(3 = — = max 

7 



/i + - /i_ p 



1 



) 



(17) 



The critical index (3 depends only on the nonlinearity index p and on the dimensionless 
velocity v = V/Vf 



This main result of our letter can be physically interpreted as follows. The exponent 
(3 is determined solely by the nonlinearity index p if the throughflow velocity is much 
larger than the front velocity (v large). Here the field in the plateau domain (see Fig. [1]) 
is effectively uncoupled from the source, and the saturation of the instability is due to 
the local nonlinearity at the source. For a small throughflow velocity (v close to one) 
the plateau state interacts with the source via the tail. Due to this "remote control," 
the field at the source is saturated more efficiently than due to nonlinearity, here the 
exponent (3 is determined solely by the form of the tail, which depends on the velocities 
ratio v. 

Below we check formula (ITS]) with direct numerical simulations of model (JTJ). A 
stationary vortex ([12]) with u = and R = 1 was imposed on a constant flow with 
V — 1. Keeping the diffusion constant fixed Dq = 0.3, for different field growth rates 
a we have found, from the linearized equations, the critical vortex intensities W cr at 
which the global mode becomes first unstable. Then we solved full nonlinear equations 
close to criticality and found the exponent (3 according to (ITT]) . The stationary problem 
was solved with a finite difference method in a domain < x < 60, < y < 40 with 
periodic boundary conditions in y and conditions u(0) = 0, f^(60) = 0. The results 
are presented in Fig. [BJ they are in very good agreement with the theoretical prediction 
( TTTj) . (|T5]) . Figure [7J shows the example of the stationary mode appearing beyond the 
instability threshold. A similar analysis performed for the quasi-one-dimensional flow 
( Hip also provides very good agreement with formula (1181) . 




(18) 
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Figure 7. The active field behind a vortex with W = 4 placed at x = y = 0, for 
V = 1, a = 0.5, D = 0.3, p = 2. 

4. Conclusions 

We have described the mixing-induced transition from a convectively unstable active 
field in an open flow to a persistent global mode. Our theoretical approach bases on the 
representation of the mixing region as a domain with locally enhanced effective diffusion. 
For a nonlinear regime close to criticality we have derived the critical exponent (3 (ITBl) 
that depends only on two parameters of the system: the dimensionless flow velocity 
v normalized by that of the front, and the nonlinearity index p. For large velocities 
the critical exponent depends only on the system's nonlinearity, which means a local 
in space saturation of the instability. For small velocities the exponent is a function of 
velocity, here the growing downstream tail of the active field imposes the saturation. 
Notably, this prediction of the one-dimensional theory is in a good accordance with two- 
dimensional calculations. The obtained results are independent of the geometry and the 
nature of the mixing region and for these reasons are expected in different systems and 
at different scales. The generic nature of our findings indicates that turbulent mixing 
can play a key role in open flows that involve active chemical and biological processes. 
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